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• Background and Aims Current understanding of stomatal development in Arabidopsis thaliana is based on 
mutations producing aberrant, often lethal phenotypes. The aim was to discover if naturally occurring viable phe- 
notypes would be useful for studying stomatal development in a species that enables further molecular analysis. 

• Methods Natural variation in stomatal abundance of A. thaliana was explored in two collections comprising 
62 wild accessions by surveying adaxial epidermal cell-type proportion (stomatal index) and density (stomatal 
and pavement cell density) traits in cotyledons and first leaves. Organ size variation was studied in a subset 
of accessions. For all traits, maternal effects derived from different laboratory environments were evaluated. 
In four selected accessions, distinct stomatal initiation processes were quantitatively analysed. 

• Key Results and Conclusions Substantial genetic variation was found for all six stomatal abundance-related 
traits, which were weakly or not affected by laboratory maternal environments. Correlation analyses revealed 
overall relationships among all traits. Within each organ, stomatal density highly correlated with the other 
traits, suggesting common genetic bases. Each trait correlated between organs, supporting supra-organ control 
of stomatal abundance. Clustering analyses identified accessions with uncommon phenotypic patterns, suggesting 
differences among genetic programmes controlling the various traits. Variation was also found in organ size, 
which negatively correlated with cell densities in both organs and with stomatal index in the cotyledon. 
Relative proportions of primary and satellite lineages varied among the accessions analysed, indicating that dis- 
tinct developmental components contribute to natural diversity in stomatal abundance. Accessions with similar 
stomatal indices showed different lineage class ratios, revealing hidden developmental phenotypes and showing 
that genetic determinants of primary and satellite lineage initiation combine in several ways. This first systematic, 
comprehensive natural variation survey for stomatal abundance in A. thaliana reveals cryptic developmental 
genetic variation, and provides relevant relationships amongst stomatal traits and extreme or uncommon acces- 
sions as resources for the genetic dissection of stomatal development. 

Key words: Natural variation, Arabidopsis thaliana, stomatal abundance, pavement cell, stomatal lineage, 
satellite stomatal lineage, epidermis, development. 



INTRODUCTION 

The potential surface available for regulated gas exchange 
between plants and the atmosphere is set by stomatal number 
and distribution in the aerial epidermis. In Arabidopsis thaliana, 
stomata differentiate gradually during organ development, 
through a series of stereotyped yet flexible cell division and 
fate acquisition events. Stomatal abundance in different plant 
surfaces and environments is regulated, resulting in variable sto- 
matal numbers and distribution patterns in mature organs 
(Bergmann and Sack, 2007; Casson and Hetherington, 2010, 
and references therein). This suggests that overlapping, partly 
redundant developmental pathways involving many genes 
must operate to produce a diversity of stomatal patterns and 
numbers while guaranteeing their functionality. Dissecting 
such gene circuits has only just begun, and several positive 
and negative regulators of stomata differentiation have been 



identified genetically and molecularly (reviewed by Bergmann 
and Sack, 2007; Nadeau, 2009; Dong and Bergmann, 2010). 

The first recognizable stomata developmental event in 
A. thaliana Col-0 (see Fig. SI in Supplementary Data, avail- 
able online) is the asymmetric division of a protodermal cell 
[meristemoid mother cell (MMC)], termed entry division, 
that initiates a stomatal cell lineage; the smaller product, the 
meristemoid (M), undergoes up to three sequential asymmetric 
amplifying divisions, oriented in an inward spiral that places 
the M in the centre of a recognizable structure made by the 
larger division products (Bergmann and Sack, 2007). While 
the central M differentiates into a guard mother cell, which 
divides symmetrically to make a stoma, each of the larger 
cells can differentiate into a pavement cell or can become a 
MMC, experience an asymmetric entry division and start a sat- 
ellite stomatal lineage (Bergmann and Sack, 2007). This is 
termed a spacing division because it puts at least one non- 
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stomatal cell between the primary (also termed planet; Lucas 
et al, 2006) and the satellite stomata. Spacing divisions 
must involve cell-cell signalling events that hinder the devel- 
opment of stomata in contact, while amplifying divisions may 
also rely on unequal distribution of stomatal fate determinants 
and, thus, be regarded as executing a lineage-based pro- 
gramme. Therefore, stomatal development is an iterative 
process involving cell lineage and cell interaction-based 
processes. 

The genetic and molecular dissection of stomatal develop- 
ment is mostly based on severe, aberrant phenotypes produced 
by induced (often loss-of-function) mutations. Such alleles, 
which impede stomata formation or produce stomatal cluster- 
ing, have identified a large suite of stomata developmental reg- 
ulators. Positive regulators are needed for stomata lineage 
initiation and development, while negative regulators are 
determinant for enforcing correct spacing (reviewed by Dong 
and Bergmann, 2010; Rowe and Bergmann, 2010). A few 
studies found some genes whose loss-of-function influence 
stomatal abundance without altering normal stomata develop- 
ment (Zhang et al, 2008; Dong and Bergmann, 2010, and 
references therein), but little is known about how diversity in 
functional, non-aberrant stomatal patterns could be set or 
how primary and satellite lineages contribute to final stomatal 
abundance. 

The quantitative traits currently used to score stomatal abun- 
dance are stomatal index (SI), which measures the proportion 
of epidermal cells that are stomata, and stomatal density 
(SD) or number of stomata per area unit. SI and SD are the 
result of cell division patterns and of cell differentiation and 
expansion during organ growth (Geisler et al, 1998; Geisler 
and Sack, 2002). SD depends on stomatal number and on the 
size and number of non-stomatal epidermal cells (mostly pave- 
ment cells), while SI depends solely on cell-type proportion, 
regardless of cell size, and therefore both traits provide comp- 
lementary information on final stomatal abundance and 
pattern. In developmental terms, the more widely examined 
organs for SI and SD are the cotyledon and first true leaf. 
Though both contribute little to total transpiration and photo- 
synthesis in the adult plant, they are crucial in the first 
stages of plant development, when transpiration is needed 
for cell expansion-mediated plant growth while seedlings 
still have a shallow root system for soil water uptake. 
Arabidopsis thaliana seed resources are mostly depleted 
3 days after germination (Penfield et al, 2005; Graham, 
2008), and seedlings rely on photosynthesis as the sole 
source of nutrients and energy for cell division and dry 
matter increase. Therefore, in early postembryonic develop- 
ment, plant survival also depends on the compromise 
between photosynthesis and transpiration, and seedling stoma- 
tal abundance and distribution are probably under selective 
pressure in natural environments. 

Naturally occurring variation among wild genotypes is an 
alternative genetic resource for studying stomatal develop- 
ment. Although many genes involved in stomatal development 
are remarkably conserved during plant evolution (Peterson 
et al., 2010), natural selection under diverse conditions may 
have resulted in natural variants that incorporate adaptive 
responses of stomata developmental gene networks to environ- 
mental cues. Indeed, natural variants in poplar and rice are 



being used to address the genetic basis of stomatal abundance 
in these taxa (Ferris et al, 2002; Laza et al, 2009). In 
A. thaliana, natural variation of complex biological processes 
has been successfully studied at the molecular level, and 
natural alleles with a broad range of quantitative effects have 
been isolated (Koornneef et al, 2004; Alonso-Blanco et al, 
2009; Lefebvre et al, 2009). Variation in the stomatal 
density response to CO2 doubling has also been reported for 
a number of A. thaliana accessions (Woodward et al, 2002). 
However, a detailed record of natural stomatal-related trait var- 
iants with traceable accessions is not yet available in this 
model species. 

MATERIALS AND METHODS 

Plant material and growth conditions 

The Arabidopsis thaliana wild genotypes studied (see Table S 1 
in Supplementary Data, available online) included 47 acces- 
sions of the Versailles 48 nested Core Collection (McKhann 
et al, 2004) and 18 accessions analysed by Clark et al. (2007) 
for sequence polymorphisms. Ler was excluded from the orig- 
inal Clark's collection because it carries non-natural alleles 
derived from a fast-neutron mutagenesis (Redei, 1962). These 
collections were built to maximize global genetic diversity in 
A. thaliana. Col-0 (N1092) was used as a laboratory reference 
accession for both collections. Three accessions were common 
to the two collections; thence, this study includes 62 different 
wild accessions plus Col-0. The INRA set seeds were obtained 
from the National Institute of Versailles' s Agronomic 
Research (INRA, France), and the remaining accessions were 
provided by the Notthingham A. thaliana Stock Centre 
(NASC). Seeds were stratified at 4 °C in darkness for 3 d, then 
sown on pots containing a 2 : 1 : 1 mixture of soil (Prohumin, 
Klasmann-Deilmann 50-50), perlite and vermiculite, and 
grown in controlled chambers (Conviron MTR30) under a 16-h 
photoperiod at 21 + 1 °C, 60-70 % relative humidity and 
150 + 20 |xmol m -2 s _1 irradiance. Data were obtained from 
four or five plants for each accession (INRA and Clark sets, 
respectively), simultaneously grown in two random blocks. 
Accessions of INRA and Clark sets were grown separately. 

Trait scoring 

The mature adaxial epidermis of cotyledons and first leaves of 
accessions were scored for stomata and pavement cell abun- 
dance. The adaxial epidermis was selected because it has 
simpler cellular patterns and fewer cell numbers than the 
abaxial epidermis. Cotyledons and first leaves of the same indi- 
viduals were analysed. To ensure that organ growth was com- 
pleted, cotyledons and first leaves were harvested 2 and 7 d 
after bolting, respectively. Sixteen accessions (see Table SI in 
Supplementary Data) did not flower in the present conditions, 
and their organs were collected 1 week later than those of the 
latest bolting accession in the assay. Surface replicas were 
obtained with dental resin (Geisler et al, 2000) and micro- 
graphed with a Leica® IRB microscope equipped with a 
Leica® DC300F camera. Images were digitally processed with 
Adobe Photoshop CS3 (Adobe Systems Inc.). Epidermal cell 
counts of each individual were an average from two 
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0-327-mm 2 areas centred along the organ apical-to-basal and the 
median-to-margin axes. Cells having at least 25 % of their 
surface inside the sampling area were scored. In leaves, this 
area excluded trichomes and their surrounding socket cells. 
Stomatal and pavement cell densities (SD and PD) were calcu- 
lated as number of stomata or pavement cells per area unit (cell 
number mm -2 ), respectively, and stomatal index (SI) as percen- 
tage of epidermal cells that were stomata. 

Organ areas were measured in the same individuals scored 
for epidermal traits. Organ replicas were photographed with 
a Leica® Stereomicroscope equipped with a Leica® 
-DC300F camera and areas determined with Leica® IM50 
image manager v. 4.0 software. 

Environmental maternal effects assay 

The maternal environment (environmental growth con- 
ditions of the mother plants of a seed progeny; ME) can 
have effects on offspring phenotypes (Donohue, 2009). Since 
ME are often non-uniform across different experiments, to 
find out if the traits under study were affected by differences 
in laboratory ME, tests were carried out. Several separate 
seed harvests were obtained for a number of accessions. 
Seven accessions representing the range of stomatal abundance 
variation in the Clark set (Bur-0, Col-0, Cvi-0, Nfa-8, 
Shakdara, Ts-1 and Van-0) were grown and selfed from the 
seed lot used for trait description, in two independent exper- 
iments (referred as ME1 and ME2) conducted in the same 
chamber at different times. Additionally, the four accessions 
present in both the Clark and the INRA collections (Bur-0, 
Col-0, Cvi-0 and Shakdara) were grown and selfed from the 
INRA collection seed stocks in a chamber different from that 
used for ME1 and ME2 (ME3). Both chambers were 
Conviron MTR30 and plants in the three ME were grown 
using the same environmental setting described above. 
Subsequently, seeds derived from a single plant for each acces- 
sion and ME were simultaneously grown, and six individuals 
per accession and ME were scored for the various traits. Data 
were analysed in two sets: one set contained four accessions 
[genotypes (G)] in three maternal environments (4G x 3ME), 
and the other had seven accessions in two ME (7G x 2ME). 

Analysis of primary and satellite stomatal lineages 

Seeds harvested from plants grown simultaneously (i.e. 
under the same ME) were grown as described above. 
Developmental homogeneity was established by monitoring 
germination every 6 h. Radicle emergence set germination 
time point at 0, and ten synchronous seedlings per accession 
were collected 48 h later. Specimens were fixed in ethanol : 
acetic acid 9 : 1 (v/v), dehydrated through ethanol : water 
series, rehydrated, mounted in Hoyer's medium (Liu and 
Meinke, 1998), and inspected with differential interference 
contrast (DIC) (Eclipse 90i microscope and Nikon 
DXM1200C camera). Epidermal cell types were recorded in 
two fields per cotyledon, avoiding the central zone and 
edges, and represented about 70 % of the organ surface. Cell 
types were identified accordingly to Zhao and Sack (1999) 
and Geisler and Sack (2002). Primary and satellite stomatal 
lineages were scored as previously defined (Berger and 



Altmann, 2000; Kutter et al, 2007; Zhang et al, 2008). No 
distinction was made between secondary and higher-order sat- 
ellite lineages. Ten additional individuals per accession were 
grown to maturity and their stomatal index determined. 

Statistical analysis 

Cell density values were log e transformed, while cell index 
values were arcsin-root transformed to improve normality of 
distributions and homogeneity of variances. None of the out- 
comes and conclusions changed when using the original 
data, and therefore most descriptions are based on untrans- 
formed data to simplify interpretation. INRA and Clark sets 
are described separately because two-way ANOVA and f-test 
analyses showed that some trait values in the four accessions 
common to both sets (Bur-0, Col-0, Cvi-0 and Shakdara) dif- 
fered significantly between assays, probably due to slight 
environmental differences resulting from their independent 
growth. For each trait, the amount of variation was calculated 
by the coefficient of variation (CV = standard deviation/ 
average) and by the fold change (maximum accession value/ 
minimum accession value). Broad-sense heritability (h 2 ) was 
calculated as: 

h 2 = V C /(V C + V E ) 

where Vq (genetic variance) is the among-genotype (acces- 
sion) variance component, and V E (environmental variance) 
is the residual error variance component estimated by 
restricted maximum likelihood (REML) analysis (Lynch and 
Walsh, 1998). Genetic correlations between traits were esti- 
mated by Pearson's and Spearman's tests using accession 
mean trait values (family mean of the selfmg offspring, 
which carry practically the same genotype; Lynch and 
Walsh, 1998). Similar relationships and conclusions were 
achieved from both analyses (unless indicated) and, hence, 
only Pearson coefficients are reported. Fisher's z-test 
(Zar, 1996) was used to identify correlations that differed sig- 
nificantly between the two accessions sets. Differences 
between mean stomatal lineage indices (arcsin-root trans- 
formed) of selected accessions were tested by Student's 
f-tests. Environmental maternal effects on traits were evaluated 
by a two-factorial analysis of variance (ANOVA), with geno- 
type (accession) and maternal environment as fixed factors. 
All statistical analyses were performed with the SPSS 11-0 
package (SPSS Inc., Chicago, IL, USA). To identify and clas- 
sify phenotypic diversity within the accession sets, hierarchical 
cluster analysis was carried out with MeV v4.2 (TM4 
Microarray Software Suite, Boston, MA, USA) using standar- 
dized data to perform an average linkage clustering based on 
uncentred Pearson's correlation as the distance metric. 

RESULTS 

Natural variation in stomatal abundance 

Stomatal abundance traits in 62 wild genotypes of A. thaliana, 
grouped in two sets (namely, the INRA and the Clark sets), 
which were grown independently and therefore described 
separately, were analysed (see Materials and methods, and 
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Table SI in Supplementary Data, available online). Col-0 was 
included in both sets as a reference genotype. Accessions were 
examined for six traits: stomatal index, stomatal density and 
pavement cell density in cotyledons (C) and first leaves (L) 
in the adaxial epidermis (namely SIC, SDC, PDC SIL, SDL 
and PDL). Pavement cell density was included to address its 
impact on stomatal density variation. Representative sampled 
epidermes for accessions with distinct stomatal abundance 
are shown in Fig. 1. No aberrant phenotypes were found in 
any of the natural accessions studied (not shown). 

All traits displayed considerable and continuous variation 
among the genotypes scored in each set (Fig. 2, and Table 
S2 in Supplementary Data). Stomatal indices showed the 
lowest variation both in cotyledons and first leaves (1-3- to 
1-7-fold). In cotyledons, SI varied 1-5-fold, ranging from the 
lowest values in Can-0 and Ts-1 to the highest values in 
Bur-0 and Van-0 in the INRA and Clark sets, respectively 
(see Table SI in Supplementary Data). Bur-0, represented in 
both sets, also had the second maximum SIC within the 
Clark group. Remarkably, Col-0 cotyledons, commonly used 
to study stomatal development, showed one of the highest sto- 
matal index values, ranking right after Bur-0 in the two data 
sets. In leaves, minimum stomatal index values occurred for 
Can-0 and Nfa-8, while maximums corresponded to 
Ishikawa (followed by Bur-0) and Rrs-7 (followed by Bur-0 
and Van-0) in the INRA and Clark sets, respectively. 




Fig. 1. Representative epidermes of A. thaliana accessions with distinct 
stomatal abundance. DIC micrographs of mature adaxial cotyledon epidermis 
(21 dpg) are shown for four accessions. Stomatal index and density values are 
low for Ts-1 (A) and Sp-0 (B), and high for Col-0 (C) and Bur-0 (D). Stomata 
have been coloured in blue. Scale bars = 100 u,m. 




The largest variation was observed in stomatal densities, 
which varied between 2-1- and 3-5-fold (see Table S2 in 
Supplementary Data). Sakata and Got-7 had the lowest SD 
values in cotyledons, and Can-0 and Tsu-1 in leaves (INRA 
and Clark sets, respectively). The highest SD corresponded 
to Jm-0 and Van-0 cotyledons, and Bur-0 and Van-0 leaves 
(INRA and Clark sets, respectively). 

Pavement cell density showed an intermediate amount of 
variation in both organs (ranging between 1-9- and 2-3-fold). 
Accessions setting the upper and lower limits of PDC were 
the same that flanked the SDC variation range (Sakata and 
Got-7 had the lowest PDC, and Jm-0 and Van-0 the highest 
one). In leaves, Ishikawa and Tsu-1 showed the lowest PDL, 
while Jm-0 and Bor-4 had the highest PDL values (INRA 
and Clark sets, respectively). 

Mean, minimum and maximum values for all traits were 
always higher in leaves than in cotyledons (Table S2 in 
Supplementary Data), suggesting a supra-organ regulation of 
stomatal abundance and pavement cell density, likely related 
to heteroblasty (Tsukaya et al, 2000). Ratios between equival- 
ent traits in leaves and cotyledons were calculated for each 
accession (Table SI in Supplementary Data). Few exceptions 
to the general behaviour were found, ratios varying from 1-3 
to 3-2 for SD, from 11 to 2-3 for PD and from 1-0 to 1-6 
for SI. Remarkably, Col-0 is one of six accessions with 
nearly identical SI in leaves and cotyledons. 

The six traits showed high to moderate broad-sense herit- 
abilities (h 2 ), which ranged from 0-33 to 0-71 (Table S2 in 
Supplementary Data). Overall, stomatal and pavement cell 
densities had similar heritability in cotyledons and leaves, 
with higher h 2 values than for stomatal index. Therefore, our 
data show that there is substantial genetic variation among 
accessions for all six traits. 

Genetic correlations among stomatal abundance traits 

To test for co-ordinated regulation of stomatal abundance in 
cotyledons and leaves and for relationships among epidermal 
traits within each organ, correlation coefficients between all 
pair combinations of the six traits were estimated (Fig. 3A-E, 
and Table S3 in Supplementary Data). Stomatal and pavement 
cell densities showed the strongest correlations, both in cotyle- 
dons and leaves (r = 0-73-0-93; P< 10 ), indicating that 
variation in pavement cell size is a major cause for SD variation 
among accessions. In addition, SI and SD positively correlated 
in both organs, with larger correlations in cotyledons (cotyle- 
dons: r=0-86, P< 10" 5 ; leaves: r=0-78, P< 10" 4 ). The 
relationship between SI and PD differed between organs; in 
cotyledons it was moderate in the two accession sets (r = 
0-49-0-61; P<5x 10" 3 ), and in leaves it was weak in the 
Clark set (r = 0-5, P < 3 x 10~ 2 ) or absent in the INRA set. 
In addition, correlations between cotyledon and leaf values for 
SI, SD and PD traits were moderately positive (Fig. 3A, and 
Table S3 in Supplementary Data). 

To identify common and distinct phenotypic patterns under- 
lying the observed correlations between trait pairs, clustering 
analyses were carried out using an uncentred Pearson correlation 
with average linkage metric (Fig. 3F, G). Dendrograms show, in 
a heat-map format, the standardized trait values (z-scores) for 
each accession, represented as a rank of divergence for the 
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Fig. 2. Frequency distributions of stomatal and pavement cell-abundance traits. Each panel shows a histogram for the 48 accessions of the NRA set (IS; upper 
part) and the 19 accessions of the Clark set (CS; lower part). Cotyledons and first leaves were scored for stomatal index [SIC (A) and SIL (D), respectively], 
stomatal density [SDC (B) and SDL (E), respectively] and pavement cell density [PDC (C) and PDL (F), respectively]. Col-0 mean values are indicated 

by a black arrow. 



trait mean in the set. In both collections, accessions were 
assigned to two main clusters: one with prevalently 
low-to-medium values (cluster type I) and another containing 
accessions with mostly high-to-medium values (cluster type 
II). Both cluster types included a group with all trait values 
ranking similarly, either low-moderate (type la; 36 % of acces- 
sions) or high-moderate (type Ila; 13 % of accessions). These 
two main clusters also included groups with trait values differ- 
ing among organs, either standing higher (lb and lib; 6 % and 
10 %, respectively) or lower (Ic; 4 %) in cotyledons than in 
leaves. In the INRA collection, cluster II included two 
additional groups with lower (lie; 19 % of accessions) or 
higher (Hd; 10 % of accessions) SI values than other traits. 
Therefore, about 49 % of accessions showed a phenotypic 
pattern with strong relationships among all trait values (la and 
Ila), in agreement with the correlations found among traits 
(Fig. 3A). The remaining groups contained phenotypes that 



deviated from these general correlations and showed, for 
instance, leaf or cotyledon stomatal index values differing 
from cell density values (see Discussion). Thus, clustering ana- 
lyses revealed natural prevalent phenotypic patterns, and also 
accessions with distinct uncommon phenotypes. 

Relationship between organ size and stomatal abundance 

Organ size regulation is expected to integrate mechanisms 
controlling stomatal lineage initiation and proliferation. To 
explore relationships between organ size and stomatal abun- 
dance, cotyledon and first leaf areas were measured in the 
INRA set. Correlations between cell and organ traits were esti- 
mated (see Fig. 4, and Tables SI, S4 and S5 in Supplementary 
Data). Substantial variation was found for cotyledon and leaf 
area (Fig. 4A, B). Notably, Sakata displayed extremely large 
organs, appearing as an outlier of the observed continuous 
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Fig. 3. Correlation analysis among epidermal cell-abundance traits. (A) Diagram showing selected correlations between traits in cotyledons and first leaves. 
Double-headed arrows connect significantly related traits. Pearson's coefficient and significance level (**, P< 0-001; *, P<001) are indicated for the 
INRA set (bold) and the Clark set (normal font). (B-E) Scatter plots showing the relationship between stomatal density and stomatal index or pavement cell 
density, in cotyledon (B, D) and first leaf (C, E). Closed circles and open circles correspond to INRA and Clark accessions, respectively. The Col-0 positions 
are indicated with red circles. (F, G) Heat-map visualization of a hierarchical cluster analysis of the cell-type abundance phenotypes at INRA (F) and Clark (G) 
accession sets. Each row represents one accession and each column represents one trait. Red and green rectangles indicate, for each trait and accession, mean 
values higher or lower than the set average; colour intensity scales are given in the top of the panel. Group designation is shown on the right and the corresponding 

nodes are highlighted in blue on the left (see text for details). Col-0 is highlighted in red. 



distribution. Hence, organ area data are presented with and 
without Sakata (Fig. 4C-F and Tables S4 and S5 in 
Supplementary Data). Broad-sense heritability estimates were 
high in both cases (ranging from 0-59 to 0-76; Table S4 in 
Supplementary Data), indicating that there is substantial 
genetic variation for these traits among natural accessions. 

Relationships among cellular and organ size traits were ana- 
lysed by Spearman correlation (Table S5 in Supplementary 
Data), since organ area data did not meet parametric test assump- 
tions required for Pearson correlations. Cotyledon and first leaf 
areas were positively correlated (p = 0-48 and 0-44; P < 0-05; 
with and without Sakata, respectively). Both showed negative 
correlations with their respective SD and PD (p between -0-29 
and -0-59; P < 0 05), particularly in cotyledons. Organ area 
and SI were negatively correlated in cotyledons (p= -0-3; 



P < 0 05), but not in leaves. Total pavement cell number for all 
accessions was estimated from their PD and organ area values 
(not shown). Cell number and organ area showed a strong positive 
correlation, which in leaves was higher than the relationship with 
any other trait (Table S5 in Supplementary Data), as previously 
reported (Granier et al, 2000; Cookson et al, 2005, 2007). 
Hence, in this sample of accessions, larger organs often had 
more and larger pavement cells (leading to lower stomatal den- 
sities) and vice versa, showing a relationship between organ 
and pavement cell number and size. All together, the present 
results show that variation in cell number and size contributes 
to genetic diversity in cotyledon and first leaf size. However, 
accessions with larger cotyledons tended to have a lower pro- 
portion of stomata (SI), suggesting a distinct genetic link 
between cotyledon size and stomatal lineage development. 
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Fig. 4. Natural diversity in cotyledon and first leaf area and their correlations with epidermal cell-abundance traits. (A, B) Histograms depict the frequency of 
accession mean values for the cotyledon area (CA) and the first leaf area (LA) in the INRA set. Values of the reference Col-0 and the extreme Sakata accessions 
are indicated (note that area scale is discontinuous). Inset: box plots for CA and LA where the upper and lower limits of the boxes represent the third and first 
quartiles and the bars inside the boxes indicate the medians; open circles indicate outlying accessions (values 1-5—3 times the interquartile range) and closed 
circles highlight the extreme Sakata value (above 3 times the interquartile range). (C-F) Scatter plots depicting relationships between (C) cotyledon and first 
leaf areas, and between cotyledon area and (D) adaxial stomatal density (SDC), (E) pavement cell density (PDC) and (F) stomatal index (SIC). Col-0 is high- 
lighted with circles. Inside each panel, the corresponding Spearman coefficient p and P-value are indicated. Sakata values were excluded in the scatter plots. 



Environmental maternal effects on epidermal cell-type abundance 
and organ size traits 

Two complementary data sets were obtained to evaluate 
whether small differences in laboratory maternal environment 
(ME), expected to occur among separate plant growth exper- 
iments, would affect the traits surveyed in this study. One 



set assessed effects of three MEs on four accessions [geno- 
types (G)], while the other used two MEs on seven accessions 
(4G x 3ME and 7G x 2ME, respectively; see Materials and 
methods). Each ME refers to a separate experiment where 
plants of the various accessions were grown simultaneously 
from seed germination to seed production. Plants from all 
maternal sources were simultaneously grown and the traits 
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scored. ME effects were estimated by two-factorial ANOVA, 
with genotype (G; accession) and ME as factors (Table S6 
in Supplementary Data). A few significant effects of ME 
(P < 0-05) were only detected in the 4G x 3ME data set. 
MEs had no effect on stomatal abundance traits, except for a 
marginal influence on leaf SD (R 2 = 2-9). Minor effects 
were also found on leaf PD (R" = 7-09). However, MEs and 
genotype x ME interactions show significant effects on 
organ sizes, being larger in cotyledons (ME, R 2 = 10-4; G x 
ME, R 2 = 20-3) than in leaves (ME, R 2 = 3-9; G x ME, 
R 2 = 13-24). 



Natural variation in stomata developmental pathways 

To address if natural variants differ in the contribution of dis- 
tinct developmental processes affecting stomatal index in 
mature organs, two quantitative processes were analysed: the 
proportion of primary lineages stemming from protodermal 
cells, and the proportion of satellite lineages that arise from non- 
stomatal cells in primary lineages (Fig. SI in Supplementary 
Data). The analyses were restricted to cotyledons in four acces- 
sions with extreme stomatal index values (Fig. 3F, G): Ts-1 and 
Sp-0, with very low SI, and Col-0 and Bur-0 as accessions with 
the highest SI values in the two collections. Adaxial cotyledon 
epidermes were examined at two time points: 48 h post- 
germination (hpg), the shortest time amenable to analyses, and 
21 d post-germination (dpg), when satellite lineages have 
appeared extensively. Lineage initiation was monitored by the 
primary lineage index (PLI) or proportion of primary stomata 
plus primary stoma precursors to total epidermal cells, and the 
satellite lineage index (SLI), or proportion of satellite stomata 
plus satellite stomata precursors (see Fig. 5 and Materials and 
methods). At 48 hpg (Table 1) the two high SI accessions 
Bur-0 and Col-0 had initiated a similar number of lineages 
(primary plus satellite), which was significantly higher than 
that of the low SI accessions Sp-0 and Ts-1. The total lineage 
index (TLI = PLI + SLI) at 48 hpg showed highly significant 
correlation with proportion of stomata (SI) in mature cotyledons 
at 21 dpg (TLI vs. SI: r = 0-99, P = 0 01), indicating that stoma- 
tal phenotypes observed at 48 hpg are predictive of mature organ 
phenotypes. In spite of this correlation, the low SI but not the 
high SI accessions showed a higher TLI at 48 h than SI at 21 
dpg (Table 1). This could be due to an increased pavement 
cell production through extended stomatal lineage amplification 
divisions and/or to more pavement cell symmetric divisions. 

Satellite lineage contribution to total stomatal lineages 
showed considerable variation among accessions, ranging 
from 10 % in Ts-1 to 33 % in Col-0 (Table 1). Although 
Bur-0 and Col-0 had different SLIs, they had a similar stomatal 
precursors index, which presented larger individual variability 
(as expected for a parameter very sensitive to small variations 
in temporal development; Fig. 5A). Regarding primary lineage 
index, Bur-0 and Col-0 showed the largest difference, while 
Bur-0 and Ts-1 displayed similar PLI values (P > 0 05). 
Therefore, two accessions with similar stomatal index 
showed a different primary lineage index (Col-0 and Bur-0), 
while accessions with a different stomatal index presented an 
equivalent primary lineage relative abundance (Bur-0 and 
Ts-1). Thus, primary and satellite lineage initiation frequencies 
can combine in different ways: low production of primary 
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Fig. 5. Diversity in distinct stomata developmental pathways in selected 
accessions. (A) Cell index on the adaxial epidermis of cotyledons 48 h after 
germination for primary stomata (PS) or precursors (PP) and satellite 
stomata (SS) or precursors (SP) in Ts-1, Sp-0, Col-0 and Bur-0. The data rep- 
resent average values of ten individuals. Note that the cell index scale is dis- 
continuous. (B-E) Representative DIC images of developing cotyledon 
epidermes in Ts-1 (B), Sp-0 (C), Col-0 (D) and Bur-0 (E). Satellite stomata 
are marked by yellow arrows and satellite precursors by blue arrowheads. 
Scale bars = 50 p,m. 



stomata with high production of satellite lineages (Col-0); 
high production of primary stomata and low satellite initiation 
(Ts-0); or moderate (Sp-0) or high (Bur-0) production of both 
primary and satellite lineages. These results reveal substantial 
natural variation in the initiation of both primary and satellite 
stomatal lineages, suggesting that the two processes are par- 
tially under independent genetic control. 

DISCUSSION 

Wild genotypes of A. thaliana show substantial genetic variation 
in epidermal cell-type abundance 

Analyses of quantitative traits related to cotyledon and first 
leaf stomatal abundance in 62 wild A. thaliana genotypes, 
evaluated in two sets, have revealed considerable natural 
genetic variation. Stomatal density is the most variable trait, 
while stomatal index showed the least variation. SI relates to 
cell division and differentiation, and its variation may result 
from relatively narrow limits that ensure a functional epidermis 
co-ordinated with the underlying mesophyll. The moderate 
pavement cell density variation found in the present study 
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Table 1. Differential contribution of primary and satellite lineage initiation to stomatal abundance in selected accessions 



Accession Primary lineages index 



Ts-1 17-9 + 04** 

Sp-0 17-1 + 0-3* 

Col-0 15-6 + 0-6 

Bur-0 18-8 + 0-5*** 



48 hpg 



Satellite lineages index Total lineages index 



2-1+0-3*** 20-0 + 0-5*** 

3.3 + 0-4*** 20-4 + 0-5*** 

7-9 + 0-7 23-4 + 0-4 

5-8 + 0-7 (*) 24-6 + 0-6 



21 dpg 

% Satellite lineages Stomatal index 



10-2+1-4*** 15-8 + 0-4*** 

16-1 + 1-7*** 16-2 + 0-3*** 

33-4 + 30 22-3 + 0-6 

23-3 + 2-5* 23-3 + 0-5 



The index of stomatal primary lineages, satellite lineages and total lineages (primary + satellite lineages), and the percentage of satellite lineages over the 
total stomatal lineages were determined in the adaxial epidermis of cotyledons 48 hpg in Ts-1, Sp-0, Col-0 and Bur-0. The stomatal index at maturity (21 dpg) 
in the same assay is provided. Significance levels of Mests with respect to Col-0 values are indicated. [***, P < 0 001; **, P < 0-01; *, P < 0 05; (*), 
P < 0-052], 



would mostly arise from cell-size phenotypes. Combinations 
of low SI and moderate PD variations most likely account 
for the higher SD diversity. 

The genetic variation found in this work suggests that phe- 
notypes of the genotypes examined are not deleterious under 
natural environments. In agreement, aberrant stomatal pheno- 
types were not found, indicating that incorrect patterns 
reduce plant fitness. To determine if the genetic variation 
observed may be involved in environmental adaptation, corre- 
lations were inspected among stomatal traits and geographic 
and historical climate data from the accession collection 
sites. Marginal positive correlations were detected between 
leaf SI and mean monthly precipitation from October to 
April (Spearman p = 0-29-0-35; P = 0-014-0-044), which 
suggests adaptive population differentiation to water avail- 
ability. However, the broad geographic distribution of the 
accessions studied and imprecise information on their original 
habitats may influence these findings. In addition, given the 
high plasticity of the traits studied (Casson and Gray, 2008; 
Lampard, 2010, and references therein), phenotypes may 
differ between laboratory conditions and natural environments, 
fading a putative adaptive basis in the underlying genetic 
variation. 

Maternal environments may influence offspring traits, 
especially those expressed early in the life cycle (reviewed 
by Donohue, 2009). It was found that cotyledon and first 
leaf cell-abundance traits were marginally or not affected by 
differences in the ME occurring in laboratory experiments. 
In agreement, A. thaliana stomatal density showed no ME 
effect during 15 generations grown at different CO2 concen- 
trations (Teng et al, 2009). As demonstrated for other 
species (Roach and Wulff, 1987), ME had a moderate but sig- 
nificant effect on cotyledon size, although genetic differences 
accounted for most of the phenotypic variation. Minor ME 
effects were also detected on first leaf size. Therefore, the phe- 
notypic differences found among accessions for cellular traits 
are mainly determined by genotypic variation and not by the 
maternal environments used. 

Relationship between epidermal cell abundance and organ size 
traits 

Correlation analyses uncovered an overall pattern of positive 
genetic relationships among stomatal and pavement 
cell-abundance traits within and between organs. The 



correlations we found most likely have a common genetic 
basis, because A. thaliana mutations in several stomata develop- 
mental genes simultaneously affect all or most traits examined 
here [notably GPA1 and ERECTA (Zhang et al, 2008; Nilson 
and Assmann, 2010; van Zanten et al, 2010); other genes are 
reviewed by Dong and Bergmann (2010)]. Nonetheless, such 
trait relationships may also originate from trait co-evolution 
(Armbruster and Schwaegerle, 1996). 

Relationships between all traits for each organ suggest that 
shared genetic networks control cell-type proportion and 
density at the organ level, although they appear more closely 
related in cotyledons than in leaves. These results are consist- 
ent with the proposed mechanism co-ordinating cell prolifer- 
ation and expansion at the organ level (Tsukaya, 2006; 
Fujikura et al, 2007, 2009; Tsukaya, 2008; Micol, 2009). 
Correlations were also found for each trait between organs. 
In agreement with these correlation patterns, cluster analyses 
identified two consistent groups of accessions where all six 
cell-abundance traits were either low (group la) or high 
(group Ha), indicating that similar genetic networks control 
cell-type abundance in both organs. These two groups 
account for nearly half of the accessions, suggesting a mechan- 
ism for supra-organ co-ordination of stomatal abundance in 
A. thaliana. Interestingly, some accessions deviated from the 
prevalent trait correlations. Clustering analysis singled out 
groups with uncommon phenotypic patterns of opposite SI 
and cell densities values (e.g. Ishikawa and Jm-0) or having 
weak or absent relationships for all traits between organs 
(e.g. C24 or Shakdara). These accessions that deviate from 
the general tendencies suggest that the genetic programmes 
controlling the various traits and their co-ordination also 
involve unshared factors. 

Cotyledon and leaf areas, which also showed significant 
genetic variation, negatively correlated with most 
cell-abundance traits in the 48 accessions examined (INRA 
set). In most accessions, large organs are built not only by 
more cells but also by larger ones (as inferred by their low 
pavement cell densities) and have lower stomatal densities; 
conversely, high stomatal and pavement cell densities (small 
cells) and lower cell numbers are the norm in small organs. 
In cotyledon, stomatal index and organ area were also nega- 
tively correlated. Overall, the present results show that wild 
accessions comply with the above-mentioned theories of 
co-ordinated cell behaviour in leaf growth (reviewed by 
Tsukaya, 2008), and suggest that, at least in cotyledons, 
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stomatal lineage development is partially linked to an organ 
level growth control. 

Some accessions carry allele combinations determining 
extremely low or high stomata proportions and cell densities, 
like Ts-1, Sp-0, Van-0 or Bur-0. The highest SD was recorded 
for Van-0, an erecta mutant (van Zanten et al, 2010). 
Loss-of-function ERECTA mutations in Col and Ler increase 
SD (Shpak et al, 2005) and total epidermal cell densities 
(Tisne et al, 2008) but do not affect SI in adult leaves 
(Masle et al., 2005), or lower it in the abaxial cotyledon epi- 
dermes (Shpak et al., 2005). It was found that Van-0 also 
has a very high pavement cell density and SI, suggesting 
differences in erecta phenotypes between epidermes or the 
presence of natural erecta modifiers in Van-0. Thus, the 
present analyses discriminated a loss-of-function allele of a 
gene involved in epidermal cell density determination, 
further revealing additional features associated with this gene. 

The wide natural genetic variation in stomatal abundance 
traits that are described here can be readily dissected by 
QTL mapping through available recombinant inbred line 
populations (Simon et al, 2008). Its combination with the 
complete genome sequences of a fast growing number of 
A. thaliana accessions (Weigel and Mott, 2009) will provide 
essential information on the mechanisms involved in this 
developmental process. Given the impact of stomatal abun- 
dance on transpiration and photosynthesis, these studies 
might provide new alleles of agronomic interest for genetic 
improvement of crop performance. 

Stomatal index encompasses cryptic genetic variation in distinct 
developmental processes 

Primary and satellite stomatal pathways have been exten- 
sively described (Bergmann and Sack, 2007; see Fig. SI in 
Supplementary Data). However, mechanisms regulating their 
proportions are poorly understood. The present analysis has 
found genetic variation in the relative frequency of primary 
and satellite lineages. A similar stomatal index could also be 
achieved with different relative frequencies. Therefore, 
primary to satellite stomata proportion is a distinct trait that 
reveals different phenotypes hidden under the same SI. This 
range of previously unknown genetic diversity suggests an 
even higher complexity of genetic factors involved in stomatal 
abundance control. 

Primary stomata abundance directly relates to entry divisions, 
while satellite stomata abundance will be affected directly by 
spacing divisions, and indirectly by amplifying and entry div- 
isions (Bergmann and Sack, 2007). Several genes are known 
to control these three asymmetric divisions (Dong and 
Bergmann, 2010), but very few seem to affect satellite lineages 
specifically. One is AGL16, whose microRNA-mediated regu- 
lation limits satellite lineage production with no apparent 
effect on primary stomata (Kutter et al., 2007). Also, two 
G-protein subunits, AGB1 and GPA1, function as mutually 
antagonistic modulators of satellite lineages (Zhang et al., 
2008). GPA1 also controls epidermal cell size (Nilson and 
Assmann, 2010). The satellite stomata fraction may largely 
depend on the cell-proliferation time window in Col-0 adaxial 
cotyledon epidermes (Geisler and Sack, 2002). The earlier cell 
division arrest in adaxial versus abaxial epidermes concurs 



with lower satellite stomata production. Therefore, changes in 
the proliferation-window length may alter satellite stomata 
abundance. Division rate differences would have a similar influ- 
ence, and both factors may also impact on higher-order 
satellization. 

Primary and satellite stomatal pathways are assumed to con- 
tribute to the phenotypic plasticity of A. thaliana stomatal 
development in response to internal and environmental cues 
(Bergmann and Sack, 2007; Casson and Gray, 2008; 
Lampard, 2010), but whether each pathway displays specific 
responses remains to be established. Unravelling the molecular 
genetic basis of the observed variation would shed light on a 
basic developmental programme, and may also provide tools 
to understand plant productivity under different environments. 

Conclusions 

A systematic data collection is provided on cotyledon and first 
leaf traits related to stomatal abundance for 62 wild Arabidopsis 
thaliana accessions analysed in two sets. This survey reveals 
substantial genetic variation in SI, SD and PD in both organs, 
and shows that these traits are largely unaffected by the labora- 
tory maternal environments used. SD shows the highest corre- 
lation with the other traits within each organ, strongly 
suggesting that these traits share genetic bases. Inter-organ cor- 
relations were also found, supporting the operation of 
supra-organ mechanisms for the control of cell-type abundance 
in the accessions surveyed. Clustering analyses showed that 
while half of the accessions had these strong relationships 
among all traits, accessions with uncommon phenotypic pat- 
terns, suggestive of differences among genetic programmes con- 
trolling the various traits, do also occur in nature. The 47 wild 
accessions surveyed had considerable variation in cotyledon 
and organ size, which were negatively correlated with cell den- 
sities. Notably, variation was also identified in two distinct pro- 
cesses during stomatal differentiation (primary and satellite 
lineage initiation), which can be hidden under stomatal abun- 
dance trait values, because accessions with a similar SI could 
have very different lineage class ratios. Thus, genetic determi- 
nants of primary and satellite lineage initiation can combine in 
several ways and produce a similar final stomatal abundance 
phenotype. In addition to revealing this cryptic diversity in a 
developmental process, this first systematic, comprehensive 
natural-variation survey for stomatal abundance in A. thaliana 
provides relevant relationships amongst stomatal traits at the 
organ and supra-organ levels. This survey also identifies acces- 
sions with extreme values or uncommon trait correlations, 
which are valuable tools for further understanding of stomatal 
development in natural germplasm. 

SUPPLEMENTARY DATA 

Supplementary data are available online at www.aob. oxford- 
journals. org and consist of the following. Figure SI: images 
illustrating the main events and cell types during stomatal 
development in A. thaliana. Table SI: accession values of 
traits analysed. Table S2: summary statistics of stomata and 
pavement cell-abundance traits in adaxial cotyledon and first 
leaf epidermes. Table S3: Correlation analysis of epidermal 
cell-abundance traits. Table S4: summary statistics of 



Delgado et al. — 



Natural variation in stomatal abundance of Arabidopsis thaliana 



1257 



cotyledon and first leaf areas analysed in the INRA set of 
accessions. Table S5: correlation analysis of organ area and 
cell-abundance traits in the INRA set of accessions. Table 
S6: Summary statistics of maternal environmental effects on 
epidermal cell abundance and organ size traits. 
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